Contaminación en el aire en la ciudad de Madrid

Authors

Gema Fernández-Avilés (Gema.FAviles@uclm.es)

Isidro Hidalgo (Isidro.Hidalgo@uclm.es)

Published

Invalid Date

Note

Los datos que se utilizan en esta historia están disponibles en el paquete CDR que puede instalarse con el siguiente comando:

install.packages("remotes")
remotes::install_github("cdr-book/CDR")

Son datos abiertos proporcionados por el portal Portal de datos abiertos del Ayuntamiento de Madrid. Concretamente, los facilitados por el Sistema Integral de la Calidad del Aire del Ayuntamiento de Madrid, que pone a disposición los datos de los contaminantes registrados por las estaciones de medición situadas en Madrid. Los datos son de frecuencia horaria por anualidades desde 2001 y los datos se actualizan de forma mensual.

Pasos:

  1. La importancia del contexto.
  2. La elección de un gráfico efectivo.
  3. El desorden es tu enemigo.
  4. Centra la atención de tu audiencia.
  5. Piensa como un diseñador.
  6. Cuenta la historia de tus datos. El archo narrativo.

1 La importancia del contexto

Es crucial proporcionar el contexto adecuado para que los lectores comprendan de dónde provienen los datos y por qué son relevantes.

Code
library(CDR) # Para los datos de temperatura
head(contam_mad)
             estaciones  id          id_name  longitud  latitud nom_mag nom_abv
1 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
2 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
3 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
4 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
5 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
6 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
  ud_med      fecha daily_mean         zona    tipo
1  µg/m3 2011-01-01  0.9375000 Interior M30 Tráfico
2  µg/m3 2011-01-02  1.0666667 Interior M30 Tráfico
3  µg/m3 2011-01-03  1.6666667 Interior M30 Tráfico
4  µg/m3 2011-01-04  1.1416667 Interior M30 Tráfico
5  µg/m3 2011-01-05  1.0416667 Interior M30 Tráfico
6  µg/m3 2011-01-06  0.4166667 Interior M30 Tráfico
Code
str(contam_mad)
Classes 'data.table' and 'data.frame':  521388 obs. of  12 variables:
 $ estaciones: chr  "E08: Escuelas Aguirre" "E08: Escuelas Aguirre" "E08: Escuelas Aguirre" "E08: Escuelas Aguirre" ...
 $ id        : chr  "E08" "E08" "E08" "E08" ...
 $ id_name   : chr  "Escuelas Aguirre" "Escuelas Aguirre" "Escuelas Aguirre" "Escuelas Aguirre" ...
 $ longitud  : num  -3.68 -3.68 -3.68 -3.68 -3.68 ...
 $ latitud   : num  40.4 40.4 40.4 40.4 40.4 ...
 $ nom_mag   : chr  "Benceno" "Benceno" "Benceno" "Benceno" ...
 $ nom_abv   : chr  "BEN" "BEN" "BEN" "BEN" ...
 $ ud_med    : chr  "µg/m3" "µg/m3" "µg/m3" "µg/m3" ...
 $ fecha     : Date, format: "2011-01-01" "2011-01-02" ...
 $ daily_mean: num  0.938 1.067 1.667 1.142 1.042 ...
 $ zona      : chr  "Interior M30" "Interior M30" "Interior M30" "Interior M30" ...
 $ tipo      : chr  "Tráfico" "Tráfico" "Tráfico" "Tráfico" ...
 - attr(*, "sorted")= chr [1:9] "estaciones" "id" "id_name" "longitud" ...
 - attr(*, ".internal.selfref")=<externalptr>
Code
library(tidyverse)

library(DT)
datatable(contam_mad[1:5, ])
Code
# Presentar las primeras líneas de los datos
head(contam_mad, 3) 
             estaciones  id          id_name  longitud  latitud nom_mag nom_abv
1 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
2 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
3 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno     BEN
  ud_med      fecha daily_mean         zona    tipo
1  µg/m3 2011-01-01   0.937500 Interior M30 Tráfico
2  µg/m3 2011-01-02   1.066667 Interior M30 Tráfico
3  µg/m3 2011-01-03   1.666667 Interior M30 Tráfico
Code
tail(contam_mad, 3) 
             estaciones  id     id_name  longitud  latitud             nom_mag
521386 E60: Tres Olivos E60 Tres Olivos -3.689731 40.50055 Óxidos de Nitrógeno
521387 E60: Tres Olivos E60 Tres Olivos -3.689731 40.50055 Óxidos de Nitrógeno
521388 E60: Tres Olivos E60 Tres Olivos -3.689731 40.50055 Óxidos de Nitrógeno
       nom_abv ud_med      fecha daily_mean    zona  tipo
521386     NOx  µg/m3 2022-04-28   32.91667 Noreste Fondo
521387     NOx  µg/m3 2022-04-29   23.83333 Noreste Fondo
521388     NOx  µg/m3 2022-04-30   35.04167 Noreste Fondo
Code
# Ver la dimensión de la tabla de datos
dim(contam_mad)
[1] 521388     12
Complejidad de los datos

Los conjuntos de datos de la calidad de aire son complejos y en algunos casos los datos no pueden utilizarse tal cual y pueden requerir consideraciones cuidadosas antes de llegar a cualquier conclusión. Debe prestarse atención a la existencia de subgrupos.

2 La elección de un gráfico efectivo

2.1 Análisis exploratorio de datos

Exercise 1 ¿Cuál es el día con mayor y menor concentración de NOx de todo el periodo?.

Code
contam_mad |> # Summary por grupo usando dplyr
  na.omit() |> # omitimos los NAs para el análisis
  filter(nom_abv == "NOx") |> # filtramos por NOx
  group_by(fecha) |> # agrupamos por fecha
  summarize(mad_mean = mean(daily_mean)) |> # promedio de las estaciones
  slice(which.max(mad_mean), which.min(mad_mean)) # seleccionamos el máximo y el mínimo
# A tibble: 2 × 2
  fecha      mad_mean
  <date>        <dbl>
1 2011-12-21   415.  
2 2020-05-10     6.33

El valor máximo, 415 µg/m3 de NOx, se observa el 21 de diciembre de 2011 y el valor mínimo, 6,32 µg/m3 de NOx, el 10 de mayo de 2020, en pleno estado de alarma.

Exercise 2 ¿Cómo son los datos de PM2.5 en la ciudad de Madrid?

Code
contam_mad |> # Summary por grupo usando dplyr
  na.omit() |> # omitimos los NAs para el análisis
  filter(nom_abv == "PM2.5") |> # filtramos por PM2.5
  group_by(id_name) |>
  summarize(
    min = min(daily_mean),
    q1 = quantile(daily_mean, 0.25),
    median = median(daily_mean),
    mean = mean(daily_mean),
    q3 = quantile(daily_mean, 0.75),
    max = max(daily_mean)
  )
# A tibble: 8 × 7
  id_name              min    q1 median  mean    q3   max
  <chr>              <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
1 Casa de Campo     0.583   5.33   8.08  9.28  11.8 177  
2 Castellana        0.125   6.08   8.83  9.96  12.3 109. 
3 Cuatro Caminos    1       6.54   9.5  10.7   13.4  59.2
4 Escuelas Aguirre  0.0417  7.25  10.3  11.5   14.3 108. 
5 Mendez Alvaro     0.0417  6.17   9.42 10.7   13.8  78.1
6 Plaza Elíptica    1.21    6.42  10.1  11.2   14   148. 
7 Plaza de Castilla 0.0417  6.29   8.92 10.1   12.1 198. 
8 Sanchinarro       1       4.47   7.27  8.60  10.5  68.5
Informes automáticos en R

Existen paquetes comno, skimr, DataExplorer y dlookr que generan informes automáticos con los principales descriptivos.

2.2 ¿Cómo han evolucionado la concentración de contaminantes en la ciudad de Madrid?

Exercise 3 Con las funciones del tidyverse represente la evolución de todos los contaminantes medidos por las estaciones de monitoreo de la ciudad de Madrid en el periodo (2011-200).

Code
plot_contam_mad <- contam_mad |>
  group_by(fecha, nom_mag) |>
  summarise(media_estaciones = mean(daily_mean, na.rm = TRUE)) |>
  ggplot(aes(x = fecha, y = media_estaciones)) +
  geom_line() +
  geom_smooth() +
  theme_minimal() +
  facet_wrap(~nom_mag, scales = "free_y")

plot_contam_mad

Vamos a mejorar el gráfico anterior para que sea más legible y efectivo.

Tu turno

Completa las partes del código señaladas por ______ o xxx para obtener el resultado propuesto.

Code
contam_mad |>
  group_by(semana = floor_date(fecha, unit = "_____"), nom_mag) |>
  summarise(media_estaciones = mean(_______, na.rm = TRUE)) |>
  ggplot(aes(x = semana, y = media_estaciones)) +
  geom_xxxxx(aes(color = nom_mag)) +
  geom_xxxxx(size = 0.5, color = "black") +
  scale_color_brewer(palette = "Paired") +
  labs(
    x = NULL, 
    y = "(µg/m3)", 
    ______ = "Evolución de partículas contaminantes en Madrid",
    ______ = "Concentración media semanal en las estaciones de medición",
    ______ = "Fuente: Portal de datos abiertos del Ayuntamiento de Madrid"
  ) +
  theme_xxxxx() +
  theme(legend.position = "none") +
  facet_wrap(~nom_mag, scales = "free_y")

2.3 Los datos faltantes

Antes de continuar, no nos olvidemos de los datos faltantes, a veces un importante problema en ciencia de datos… ¿Existe algún patrón en los NAs?

Exercise 4 ¿Cuántos NAs tengo en mi conjunto de datos contam_mad?

Code
sum(is.na(contam_mad))
[1] 15615

Exercise 5 ¿Puedo visualizar dónde están los datos faltantes por estación de monitoreo a lo largo del tiempo?

Code
na_table <- contam_mad |>
  mutate(isna = is.na(daily_mean)) |>
  ggplot(aes(x = fecha, y = id_name, fill = isna)) +
  geom_raster() +
  scale_fill_manual(values = c("TRUE" = "red", "FALSE" = "lightblue"))  

na_table

Exercise 6 ¿Qué hacemos con los NAs?

Opción 1: Eliminar los NAs

Code
contam_mad_clean <- contam_mad |>
  drop_na()

summary(is.na(contam_mad_clean))
 estaciones          id           id_name         longitud      
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:505773    FALSE:505773    FALSE:505773    FALSE:505773   
  latitud         nom_mag         nom_abv          ud_med       
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:505773    FALSE:505773    FALSE:505773    FALSE:505773   
   fecha         daily_mean         zona            tipo        
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:505773    FALSE:505773    FALSE:505773    FALSE:505773

Opción 2: Imputar NAs con el día anterior

Code
contam_mad_dia_antes <- contam_mad |>
  arrange(estaciones, nom_abv, fecha) |>
  fill(daily_mean)

summary(is.na(contam_mad_dia_antes))
 estaciones          id           id_name         longitud      
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:521388    FALSE:521388    FALSE:521388    FALSE:521388   
  latitud         nom_mag         nom_abv          ud_med       
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:521388    FALSE:521388    FALSE:521388    FALSE:521388   
   fecha         daily_mean         zona            tipo        
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:521388    FALSE:521388    FALSE:521388    FALSE:521388
Code
#imputar los na con la mediana de la estación

contam_mad_mediana <- contam_mad |>
  group_by(estaciones) |>
  fill(daily_mean, .direction = "updown")
Tu turno

Completa las partes del código señaladas por ______ o xxx para obtener el resultado esperado.

Code
summary(____(contam_mad_mediana))
 estaciones          id           id_name         longitud      
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:521388    FALSE:521388    FALSE:521388    FALSE:521388   
  latitud         nom_mag         nom_abv          ud_med       
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:521388    FALSE:521388    FALSE:521388    FALSE:521388   
   fecha         daily_mean         zona            tipo        
 Mode :logical   Mode :logical   Mode :logical   Mode :logical  
 FALSE:521388    FALSE:521388    FALSE:521388    FALSE:521388

3 El desorden es tu enemigo: calendario de contaminantes

Mantén el código limpio y organizado para facilitar su comprensión y mantenimiento. Veamos un ejemplo con el siguiente ejemplo. Una herramienta muy útil para tener una visión general de estos contaminantes es viendo el calendario como un heatmap: calendar heatmap

Preparamos un código para visualizar la concentración media de los contaminantes en Madrid a lo largo del tiempo en forma de calendario que pueda ser utilizado para cualquier contaminante.

Code
calendar_plot <- contam_mad |>
  group_by(fecha, nom_mag, nom_abv, ud_med) |>
  summarize(valor_promedio = mean(daily_mean, na.rm = T))

# Dates as factors
months <-
  seq.Date(
    from = as.Date("2022-01-01"),
    length.out = 12,
    by = "month"
  ) |> format("%B")
wdays <-
  seq.Date(
    from = as.Date("2022-05-30"),
    length.out = 7,
    by = "day"
  ) |> format("%A")

calendar_plot <- calendar_plot |>
  mutate(
    year = format(fecha, "%Y"),
    month = factor(format(fecha, "%B"), levels = months, labels = months),
    wday = factor(weekdays(fecha), levels = wdays, labels = wdays),
    week = as.numeric(format(fecha, "%W"))
  )

calendar_plot <- calendar_plot |>
  group_by(year, month) |>
  mutate(wmonth = 1 + week - min(week))

Exercise 7 Calendar heatmap para el NOx

Code
i_mag <- "Óxidos de Nitrógeno"    # Seleccionar el contaminante

fill_title <- calendar_plot |>
  filter(nom_mag == i_mag & year >= 2011) |>
  ungroup() |>
  distinct(paste(unique(nom_abv), unique(ud_med)))

calendar_plot |>
  filter(nom_mag == i_mag & year >= 2011) |>
  ggplot(aes(
    x = wmonth,
    y = reorder(wday, -as.numeric(wday)),
    fill = valor_promedio
  )) +
  geom_tile(colour = "white") +
  facet_grid(year ~ month) +
  scale_fill_gradient(low = "yellow", high = "red", ) +
  scale_x_continuous(breaks = 1:5, limits = c(0, 6)) +
  labs(
    x = "Semana del mes",
    y = NULL,
    title = paste0("Concentración de ", i_mag, " por día de la semana"),
    fill = fill_title,
    caption = "Fuente: Red de Vigilancia de la Calidad del Aire del Ayto. de Madrid"
  )

Tu turno

Completa las partes del código señaladas por ______ o xxx para obtener el calendario de PM10.

Code
i_mag <- "_______________"

fill_title <- calendar_plot |>
  filter(nom_mag == i_mag & year >= 2011) |>
  ungroup() |>
  distinct(paste(unique(nom_abv), unique(ud_med)))

calendar_plot |>
  filter(nom_mag == i_mag & year >= 2011) |>
  ggplot(aes(
    x = wmonth,
    y = reorder(wday, -as.numeric(wday)),
    fill = valor_promedio
  )) +
  geom_tile(colour = "white") +
  facet_grid(year ~ month) +
  scale_fill_gradient(low = "________", high = "red", ) +
  scale_x_continuous(breaks = 1:5, limits = c(0, 6)) +
  labs(
    x = "Semana del mes",
    y = NULL,
    title = paste0("Concentración de ", i_mag, " por día de la semana"),
    fill = fill_title,
    caption = "Fuente: Red de Vigilancia de la Calidad del Aire del Ayto. de Madrid"
  )

4 Centra la atención de tu audiencia

Exercise 8 Nos centramos en los dos contaminantes más problemáticos en la ciudad de Madrid (PM10 y NOx)

Code
contam_mad_pm10_nox <- contam_mad |>
  filter(nom_abv %in% c("PM10", "NOx"))
Code
# echo: false

plot_pm10_nox <- contam_mad_pm10_nox |>
  group_by(fecha, nom_mag) |>
  summarise(media_estaciones = mean(daily_mean, na.rm = TRUE)) |>
  ggplot(aes(x = fecha, y = media_estaciones)) +
  geom_line(aes(color = nom_mag)) +
  labs(
    x = NULL, y = "(µg/m3)", title = "Evolución semanal de partículas contaminantes (PM10 y NOx) en Madrid",
    subtitle = "Concentración media semanal en las estaciones de medición",
    caption = "Fuente: Portal de datos abiertos del Ayuntamiento de Madrid"
  ) +
  theme_bw() +
  theme(legend.position = "none") +
  facet_wrap(~nom_mag, scales = "free_y", ncol = 1)

plot_pm10_nox

Tu turno

Completa las partes del código señaladas por ______ o xxx para obtener el resultado deseado.

Code
plot_pm10_nox <- ____________  |>
  group_by(fecha, nom_mag) |>
  summarise(media_estaciones = mean(daily_mean, na.rm = TRUE)) |>
  ggplot(aes(x = ___________, y = ____________)) +
  geom_line(aes(color = nom_mag)) +
  labs(
    x = NULL, y = "(µg/m3)", 
    title = "Evolución semanal de PM10 y NOx en Madrid",
    subtitle = "Concentración media semanal en las estaciones de medición",
    caption = "Fuente: Portal de datos abiertos del Ayuntamiento de Madrid"
  ) +
  theme_bw() +
  theme(legend.position = "none") +
  facet_wrap(~nom_mag, scales = "free_y", ncol = 1)

plot_pm10_nox

Exercise 9 ¿Por qué no hacer el anterior gráfico interactivo?

La función ggplotly() de la librería plotly permite hacer fácilmente gráficos interactivos.

Code
plotly::ggplotly(plot_pm10_nox)

Exercise 10 Algunas relaciones bi-variantes y n-variantes interesantes que se pueden establecer entre los contaminantes PM10 y NOx en el periodo del estado de alarma por la pandemia de COVID-19.

Code
pm10_nox_mad_alarma <- contam_mad |>
  na.omit() |>
  filter(nom_abv %in% c("PM10", "NOx")) |>
  # período del estado de alarma
  filter(between(fecha, left = as.Date("2020-03-14"), right = as.Date("2020-06-30"))) |>
  select(estaciones, zona, tipo, nom_abv, daily_mean, fecha) |>
  pivot_wider(names_from = "nom_abv", values_from = "daily_mean", values_fn = mean)

pm10_nox_mad_alarma |>
  ggplot(
    aes(x = PM10, y = NOx, colour = tipo, size = zona)
  ) +
  geom_point()

Code
pm10_nox_mad_alarma |>
  drop_na() |>
  ggplot(aes(y = NOx, x = PM10, colour = tipo, shape = zona)) +
  geom_point() +
  geom_smooth() +
  facet_wrap(vars(estaciones))

5 Piensa como un diseñador

Antes de avanzar vamos a poner a prueba lo aprendido hasta aquí. Para ello nos centraremos en el NOx.

Code
nox_madrid <- contam_mad |>
  na.omit() |>
  filter(nom_abv == "NOx") |>
  filter(between(fecha, left = as.Date("2022-03-01"), right = as.Date("2022-03-31")))
Tu turno

Completa las partes del código señaladas por ______ o xxx para obtener el resultado deseado.

Code
nox_madrid |>
  _____(aes(x= zona, y=daily_mean)) +
  geom_xxxxx() +
  geom_xxxxx(height = 0, width = 0.01) +
  aes(x = zona, y = daily_mean, fill = zona) +
  labs(
    title = "Distribución de NOx por zona",
    x = "Zona",
    y = "Concentración de NOx (µg/m3)"
  ) +
  theme_xxxxx()

Ahora con algunas estadísticas interesantes:

Tu turno

Completa las partes del código señaladas por ______ o xxx para obtener el resultado deseado.

Code
library(ggstatsplot)

ggbetweenstats(
  data = ______,
  x = _______,
  y = _______
)

Ahora, utilizando la función geom_density_ridges() de la librería ggridges vamos a hacer un gráfico de densidad (Ridgeline) de NOx por zona.

Code
library(ggridges)
nox_madrid |>
  ggplot(aes(x = daily_mean, y = zona, fill = zona)) +
  geom_density_ridges() +
  theme_ridges() +
  labs(
    title = "Distribución de NOx por zona de calidad del aire",
    x = "Concentración de NOx (µg/m3)",
    y = "Zona"
  ) +
  theme_minimal()

Exercise 11 ¿Cómo es la relación bivariante de estos dos contaminantes?

Compruebe con un Análisis de la Varianza si los niveles de concentración de NOx dependen de las variables tipo y zona

Code
contam_mad_anova <- contam_mad |>
  filter(nom_abv == "NOx") |>
  drop_na()

anova <- aov(data = contam_mad_anova, daily_mean ~ zona + tipo)
summary(anova)
               Df    Sum Sq Mean Sq F value Pr(>F)    
zona            4  22508217 5627054  1540.1 <2e-16 ***
tipo            2   5489043 2744521   751.2 <2e-16 ***
Residuals   94975 347009104    3654                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
ggplot(contam_mad_anova, aes(x = zona, y = daily_mean, fill = tipo)) +
  geom_boxplot() +
  labs(title = "NOx por Zona y Tipo",
       x = "Zona",
       y = "Daily Mean",
       fill = "Tipo") +
  theme_minimal()

¿Qué paso la semana de la calima de marzo de 2022 con el PM en la ciudad de Madrid?

Code
particulas <- c("Partículas < 2.5 µm", "Partículas < 10 µm")

calima <- contam_mad |>
  filter(nom_mag %in% particulas &
    fecha %in% seq.Date(as.Date("2022-03-01"), by = "day", length.out = 31)) |>
  group_by(fecha, id, id_name, nom_mag, nom_abv, ud_med) |>
  summarize(valor_promedio = mean(daily_mean, na.rm = T))


max_2.5 <- calima |>
  ungroup() |>
  filter(nom_mag == particulas[1]) |>
  slice(which.max(valor_promedio))
max_10 <- calima |>
  ungroup() |>
  filter(nom_mag == particulas[2]) |>
  slice(which.max(valor_promedio))

library(ggrepel)
calima |>
  ggplot(aes(fecha, valor_promedio, colour = nom_mag)) +
  geom_jitter() +
  geom_smooth(
    method = "loess",
    span = .5,
    se = FALSE,
    show.legend = FALSE
  ) +
  scale_x_date(
    breaks = seq.Date(as.Date("2022-03-01"), by = "week", length.out = 5),
    date_labels = "%d-%b"
  ) +
  scale_color_manual(values = c("#261606", "#DD9C4A")) +
  geom_label_repel(
    data = max_10,
    mapping = aes(fecha, valor_promedio, label = paste(id_name)),
    show.legend = FALSE
  ) +
  geom_label_repel(
    data = max_2.5,
    mapping = aes(fecha, valor_promedio, label = paste(id_name)),
    show.legend = FALSE
  ) +
  labs(
    title = "Registro de partículas durante el mes de marzo 2022",
    subtitle = "Madrid",
    x = NULL,
    y = unique(calima$ud_med),
    color = NULL,
    caption = "Fuente: Red de Vigilancia de la Calidad del Aire del Ayto. de Madrid"
  ) +
  theme_minimal()

6 Cuenta la historia de tus datos. Predicción espacial del PM10 (calima del 13-17 marzo)

Code
# idw
library(mapSpain)
library(sf)

mad_sf <- esp_get_munic(munic = "^Madrid$", epsg = 4326)

marzo_pm10 <- contam_mad |>
  filter(nom_abv == "PM10" & fecha >= as.Date("2022-03-13") & fecha <= as.Date("2022-03-17")) |>
  drop_na()

madrid_estaciones_sf <- st_as_sf(marzo_pm10,
  coords = c("longitud", "latitud"),
  crs = 4326
)

library(gganimate)
ggplot(madrid_estaciones_sf) +
  geom_sf(
    data = mad_sf,
    fill = "#DD9C4A",
  ) +
  geom_sf(aes(fill = daily_mean),
    shape = 21,
    size = 5,
    alpha = .7
  ) +
  labs(fill = "PM10") +
  scale_fill_viridis_c() +
  theme_void() +
  labs(
    title = "PM10: {current_frame}"
  ) +
  transition_manual(fecha) +
  ease_aes("linear") +
  theme(
    plot.title = element_text(
      size = 12,
      face = "bold"
    ),
    plot.subtitle = element_text(
      size = 8,
      face = "italic"
    )
  )